function [b,vc] = vciqr_oid(bhat,y,d,x,z,tau)
% Compute standard errors for IQR when the equation is just identified
% using the method proposed in Powell (1986).  This implementation uses
% a Gaussian kernel and a simple rule of thumb bandwidth.
% [b,vc,J] = vciqr(bhat,y,d,x,z,tau)
% Inputs -
%   bhat - estimated coefficients from IQR
%   y - dependent variable
%   d - RHS endogenous variable
%   x - covariates (If there are no covariates, pass x = [].)
%   z - instruments
%   tau - quantile index
% Outputs -
%   b - estimated coefficients with standard errors
%   vc - covariance matrix of b
%   J - matrix in asymptotic variance formula (J'SJ) used in process testing


n = size(y,1);	      % Number of observations
x = [x,ones(n,1)];    % Add constant term
k = size([d,x],2);    % Number of regressors

vc = zeros(k,k);
b = zeros(k,2);
S = (tau-tau^2)*(1/n)*[z,x]'*[z,x]; 

e = y - [d,x]*bhat; % Generate residuals

h  = 1.364*((2*sqrt(pi))^(-1/5))*sqrt(var(e))*(n^(-1/5));            %Calculate bandwidth using Silverman's rule of thumb
Ja = (1/(n*h))*((normpdf(e/h)*ones(1,size(d,2))).*d)'*[z,x];
H  = (1/(n*h))*((normpdf(e/h)*ones(1,size([z,x],2))).*[z,x])'*[z,x];
JT = inv(H);
Jg = JT(1:size(z,2),:);
Jb = JT(size(z,2)+1:size(JT,1),:);

K  = pinv(Ja*Jg'*pinv(Jg*S*Jg')*Jg*Ja')*(Ja*Jg'*pinv(Jg*S*Jg')*Jg);
L  = Jb - Jb*Ja'*K;

M = [K ; L];

vc = (1/n)*M*S*M';
b(:,1) = bhat;
b(:,2) = (sqrt(diag(vc)));
